ESP - Geophysical Fluid Dynamics

Villanueva, Lloyd

# Needed libraries and functions
from IPython.display import HTML, display
from matplotlib.animation import FuncAnimation
import numpy as np
import matplotlib.pyplot as plt

Exercise II: Kinematics

1. A two-dimensional steady flow has velocity components u = y and v = x. Show that the streamlines are rectangular hyperbolas

\begin{equation} x^2 - y^2 = const. \end{equation}

Sketch the flow pattern and convince yourself that it represents an irrotational flow in a $90 [^\circ]$ corner.

To show that the streamlines are a rectangular hyperbolas, we will start from the streamline equation \begin{align} \frac{\mathrm{d}y}{v} & = \frac{\mathrm{d}x}{u} \\ \frac{\mathrm{d}y}{x} & = \frac{\mathrm{d}x}{y} \\ y \mathrm{d}y & = x \mathrm{d}x \\ \int y \, \mathrm{d}y & = \int x \, \mathrm{d}x \quad \quad \text{integrate both sides of the equation}\\ \frac{y^2}{2} + C_1 & = \frac{x^2}{2} + C_2 \\ y^2 & = x^2 + 2(C_2 - C1) \\ x^2 - y^2 & = -2(C_2 - C1) \quad \quad \text{combine all the constant and rearrange the equation} \\ x^2 - y^2 & = const. \quad \quad \text{equation of the rectangular parabolas} \end{align}

import numpy as np
import matplotlib.pyplot as plt

# Define the velocity components
def velocity_components(x, y):
    u = y
    v = x
    return u, v

# Generate grid points
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)

# Calculate velocity components at each grid point
U, V = velocity_components(X, Y)

# Plot streamlines
plt.figure(figsize=(8, 8))
plt.streamplot(X, Y, U, V, density=1, linewidth=1, color='blue', arrowstyle='->', arrowsize=1.5)

# Set axis labels and title
plt.xlabel('x')
plt.ylabel('y')
plt.title('Streamlines of $u = y$ and $v = x$')

# Show the plot
plt.grid(True)
plt.axhline(0, color='black',linewidth=0.5)
plt.axvline(0, color='black',linewidth=0.5)
plt.gca().set_aspect('equal', adjustable='box')
plt.show()

To show its irrotational, below equation must be satisfy \begin{equation} \nabla \times V = 0. \end{equation} For this case, $V = (u, v) = (y, x)$ and its curl is \begin{align} \nabla \times V & = 0 \\ \left( \frac{\delta}{\delta x} v - \frac{\delta}{\delta y} u \right) \hat{k} & = 0 \\ \frac{\delta}{\delta x} x - \frac{\delta}{\delta y} y & = 0 \\ 1 - 1 & = 0 \\ 0 & = 0 \end{align} Thus, the flow pattern is irrotational.

2. The velocity components in an unsteady plane flow are given by

\begin{equation} u = x / (1 + t) \end{equation}

and

\begin{equation} v = 2y / (2 + t) \end{equation}

Describe the path lines and the streamlines. Note that the path lines are found by following the motion of each particle, that is, by solving the differential equations $\mathrm{d}x/\mathrm{d}t = u(x, t)$ and $\mathrm{d}y/\mathrm{d}t = v(x, t)$, subject to $x = x_0$ at $t = 0$.

Solution:

Path lines

Solving the two differential for the path lines \begin{equation} \frac{\mathrm{d}x}{\mathrm{d}t} = u(x,t) \quad\mathrm{and} \quad \frac{\mathrm{d}y}{\mathrm{d}t} = v(x,t). \end{equation} First, solving the first differential equation, we will have \begin{align} \frac{\mathrm{d}x}{\mathrm{d}t} & = u(x,t) \\ \frac{\mathrm{d}x}{\mathrm{d}t} & = \frac{x}{1+t} \\ \frac{\mathrm{d}x}{x} & = \frac{\mathrm{d}t}{1+t} \\ \int_{x' = x_0}^{x' = x} \frac{\mathrm{d}x'}{x'} & =\int_{t'=0}^{t' = t} \frac{\mathrm{d}t'}{1+t'} \\ \ln{x} - \ln{x_0} & = \ln{(1+t)} - \ln{1} \\ \ln{x} & = \ln{(1+t)} + \ln{x_0} \\ \ln{x} & = \ln{(1+t)(x_0)} \end{align}

\begin{equation} x(t) = (1+t)(x_0) \end{equation}.

Solving for the second differential equation, \begin{align} \frac{\mathrm{d}y}{\mathrm{d}t} & = v(y,t) \\ \frac{\mathrm{d}y}{\mathrm{d}t} & = \frac{2y}{2+t} \\ \frac{1}{2} \frac{\mathrm{d}y}{y} & = \frac{\mathrm{d}t}{2+t} \\ \frac{1}{2} \int_{y' = y_0}^{y' = y} \frac{\mathrm{d}y'}{y'} & = \int_{t'=0}^{t' = t} \frac{\mathrm{d}t'}{2+t'} \\ \frac{1}{2} \ln{y} - \frac{1}{2} \ln{y_0} & = \ln{(2+t)} - \ln{2} \\ \ln{y} & = 2\ln{(2+t)} -2 \ln{2} + \ln{y_0} \\ \ln{y} & = \ln{\frac{(2+t)^{2}y_0}{4}} \end{align}

\begin{equation} y = \frac{(2+t)^2y_0}{4} \end{equation}

Next, to be able to plot the path lines, we need to find the expression of $t$ from any of the two solutions and substitute it to the other solutions \begin{align} x & = (1+t)(x_0) \\ x & = x_0 + x_0t \\ x - x_0 & = x_0t \\ t & = \frac{x-x_0}{x_0} \\ t & = \frac{x}{x_0} - 1 \end{align} Substituting it to the other equation, we will have \begin{align} y & = \frac{(2+t)^2y_0}{4} \\ y & = \frac{(2+\frac{x}{x_0} - 1)^2 y_0}{4} \\ y & = \frac{(1+\frac{x}{x_0})^2 \, y_0}{4} \end{align}

\begin{equation} y = \frac{y_0}{4} \, \left(1+\frac{x}{x_0} \right)^2 \end{equation}

The pathlines obtained is a parabola, which opening is dependent on the sign of $y_0$ (going up if (+) and going down if (-)).

Streamlines

For the streamlines, we will start from the streamline equation, \begin{align} \frac{\mathrm{d}y}{v} & = \frac{\mathrm{d}x}{u} \\ \frac{\mathrm{d}y}{2y/(2+t)} & = \frac{\mathrm{d}x}{x/(1+t)} \\ \frac{\mathrm{d}y}{2y} & = \frac{1+t}{2+t} \frac{\mathrm{d}x}{x} \\ \frac{1}{2} \int_{y' = y_0} ^{y'=y}\frac{1}{y'} \, \mathrm{d}y' & = \frac{1+t}{2+t} \, \int_{x' = x_0} ^{x'=x} \frac{1}{x'} \, \mathrm{d}x' \quad \quad \text{integrate both sides of the equation}\\ \frac{1}{2} \ln{y} - \frac{1}{2} \ln{y_0} & = \frac{1+t}{2+t} \, (\ln{x} - \ln{x_0}) \\ \ln{\frac{y}{y_0}} & = \frac{2(1+t)}{2+t} \, \ln{\frac{x}{x_0}} \\ \ln{\frac{y}{y_0}} & = \ln{ \left( \frac{x}{x_0} \right)^{\frac{2(1+t)}{2+t}}} \\ \frac{y}{y_0} & = \left( \frac{x}{x_0} \right)^{\frac{2(1+t)}{2+t}} \end{align}

\begin{equation} y= y_0 \left( \frac{x}{x_0} \right)^{\frac{2(1+t)}{2+t}} \end{equation}

For the streamlines, at $t=0$ the streamline is a just a straight line going outward from the origin which is shown the plot below. As the time progress, the degree of $x$ increases and this is observed in the simulation.

import numpy as np
import matplotlib.pyplot as plt

# Define the velocity components
def velocity_components(x, y):
    u = x
    v = y
    return u, v

# Generate grid points
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)

# Calculate velocity components at each grid point
U, V = velocity_components(X, Y)

# Plot streamlines
plt.figure(figsize=(8, 8))
plt.streamplot(X, Y, U, V, density=1, linewidth=1, color='black', arrowstyle='->', arrowsize=1.5)

# Set axis labels and title
plt.xlabel('x')
plt.ylabel('y')
plt.title('Streamlines of u = x/(1+t) and v = 2y/(2+t) at t=0')

# Show the plot
plt.grid(True)
plt.axhline(0, color='black',linewidth=0.5)
plt.axvline(0, color='black',linewidth=0.5)
plt.gca().set_aspect('equal', adjustable='box')
plt.show()

def your_u_function(x, y, t):
    u = x / (1+t)
    return u

def your_v_function(x, y, t):
    v = (2*y) / (2+t)
    return v

# Define your 2D vector field function (u, v) and time array
def vector_field(x, y, t):
    # Replace with your actual u and v functions
    u = your_u_function(x, y, t)
    v = your_v_function(x, y, t)
    return u, v

# Set up the figure and axes
fig, ax = plt.subplots(figsize=(8, 6))
ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)

# Create a grid of points
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)

# Function to update the plot for each frame
def update(frame):
    ax.clear()
    ax.set_xlim(-2, 2)
    ax.set_ylim(-2, 2)
    
    # Evaluate the vector field at each point on the grid
    U, V = vector_field(X, Y, frame)
    
    # Plot streamlines
    stream = ax.streamplot(X, Y, U, V, color='black')
    
    # Add a title or other annotations if needed
    ax.set_title(f'Time: {frame}')
    
# Create an animation
num_frames = 15  # Adjust the number of frames as needed
animation = FuncAnimation(fig, update, frames=num_frames, interval=100)

# Save the animation as a GIF
animation.save('streamlines_animation.gif', writer='pillow')

# Display the animation in the Jupyter Notebook
html_output = HTML(animation.to_jshtml())
display(html_output)

Once Loop Reflect

3. Consider the simple plane flow given by

\begin{align} u_1 & = x_1 / (1 + t) \\ u_2 & = x_2 \\ u_3 & = 0 \end{align}

Compute the velocity gradient matrix, strain and rotation tensors. Compute divergence and vorticity of the given field.

The velocity gradient matrix is defined as \begin{equation} \frac{\partial u_i}{\partial u_j} = \frac{1}{2} \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) + \frac{1}{2} \left( \frac{\partial u_i}{\partial x_j} -\frac{\partial u_j}{\partial x_i} \right) = \epsilon_{ij} + \frac{1}{2} r_{ij} \end{equation} where \begin{equation} \epsilon_{ij} = \frac{1}{2} \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) \end{equation} is the strain rate tensor and \begin{equation} r_{ij} = \left( \frac{\partial u_i}{\partial x_j} - \frac{\partial u_j}{\partial x_i} \right) \end{equation} is the rotation tensor.

Strain rate

First, we will calculate the strain rate tensor of the given flow \begin{align} \epsilon_{11} & = \frac{1}{2} \left( \frac{\partial u_1}{\partial x_1} + \frac{\partial u_1}{\partial x_1} \right) = \frac{1}{2} \left( \frac{\partial}{\partial x_1} \frac{x_1}{1+t} + \frac{\partial}{\partial x_1} \frac{x_1}{1+t} \right) = \frac{1}{2} \left( \frac{1}{1+t} + \frac{1}{1+t} \right) = \frac{1}{1+t} \\ \epsilon_{12} & = \frac{1}{2} \left( \frac{\partial u_1}{\partial x_2} + \frac{\partial u_2}{\partial x_1} \right) = \frac{1}{2} \left( \frac{\partial}{\partial x_2} \frac{x_1}{1+t} + \frac{\partial}{\partial x_1} x_2 \right) = 0 \\ \epsilon_{13} & = \frac{1}{2} \left( \frac{\partial u_1}{\partial x_3} + \frac{\partial u_3}{\partial x_1} \right) =\frac{1}{2} \left( \frac{\partial}{\partial x_3} \frac{x_1}{1+t} + 0\right) = 0 \\ \epsilon_{21} & = \frac{1}{2} \left( \frac{\partial u_2}{\partial x_1} + \frac{\partial u_1}{\partial x_2} \right) = \frac{1}{2} \left( \frac{\partial}{\partial x_1} x_2 + \frac{\partial}{\partial x_2} \frac{x_1}{1+t} \right) = 0 \\ \epsilon_{22} & = \frac{1}{2} \left( \frac{\partial u_2}{\partial x_2} + \frac{\partial u_2}{\partial x_2} \right) = \frac{1}{2} \left( \frac{\partial}{\partial x_2} x_2 + \frac{\partial}{\partial x_2} x_2 \right) = \frac{1}{2} \left( 1 + 1 \right) = 1 \\ \epsilon_{23} & = \frac{1}{2} \left( \frac{\partial u_2}{\partial x_3} + \frac{\partial u_3}{\partial x2} \right) =\frac{1}{2} \left( \frac{\partial}{\partial x_3} x_2 + 0\right) = 0 \\ \epsilon_{31} & = \frac{1}{2} \left( \frac{\partial u_3}{\partial x_1} + \frac{\partial u_1}{\partial x_3} \right) = \frac{1}{2} \left( 0 + \frac{\partial}{\partial x_2} \frac{x_1}{1+t} \right) = 0 \\ \epsilon_{32} & = \frac{1}{2} \left( \frac{\partial u_3}{\partial x_2} + \frac{\partial u_2}{\partial x_3} \right) = \frac{1}{2} \left( 0 +\frac{\partial}{\partial x_3} x_2 \right) = 0 \\ \epsilon_{33} & = \frac{1}{2} \left( \frac{\partial u_3}{\partial x_3} + \frac{\partial u_3}{\partial x3} \right) =\frac{1}{2} \left( 0 + 0 \right) = 0, \end{align}

thus, the strain rate tensor is calculated to be \begin{equation} \epsilon_{ij} = \begin{bmatrix} \frac{1}{1+t} & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{bmatrix} \end{equation}

Rotation

Now, we calculate the rotation tensor \begin{align} r_{11} & = \left( \frac{\partial u_1}{\partial x_1} - \frac{\partial u_1}{\partial x_1} \right) = \left( \frac{\partial}{\partial x_1} \frac{x_1}{1+t} - \frac{\partial}{\partial x_1} \frac{x_1}{1+t} \right) = \left( \frac{1}{1+t} - \frac{1}{1+t} \right) = 0 \\ r_{12} & = \left( \frac{\partial u_1}{\partial x_2} - \frac{\partial u_2}{\partial x_1} \right) = \left( \frac{\partial}{\partial x_2} \frac{x_1}{1+t} - \frac{\partial}{\partial x_1} x_2 \right) = 0 \\ r_{13} & = \left( \frac{\partial u_1}{\partial x_3} - \frac{\partial u_3}{\partial x_1} \right) = \left( \frac{\partial}{\partial x_3} \frac{x_1}{1+t} - 0\right) = 0 \\ r_{21} & = \left( \frac{\partial u_2}{\partial x_1} - \frac{\partial u_1}{\partial x_2} \right) = \left( \frac{\partial}{\partial x_1} x_2 - \frac{\partial}{\partial x_2} \frac{x_1}{1+t} \right) = 0 \\ r_{22} & = \left( \frac{\partial u_2}{\partial x_2} - \frac{\partial u_2}{\partial x_2} \right) = \left( \frac{\partial}{\partial x_2} x_2 - \frac{\partial}{\partial x_2} x_2 \right) = \left( 1 - 1 \right) = 0 \\ r_{23} & = \left( \frac{\partial u_2}{\partial x_3} - \frac{\partial u_3}{\partial x_2} \right) = \left( \frac{\partial}{\partial x_3} x_2 - 0\right) = 0 \\ r_{31} & = \left( \frac{\partial u_3}{\partial x_1} - \frac{\partial u_1}{\partial x_3} \right) = \left( 0 - \frac{\partial}{\partial x_2} \frac{x_1}{1+t} \right) = 0 \\ r_{32} & = \left( \frac{\partial u_3}{\partial x_2} - \frac{\partial u_2}{\partial x_3} \right) = \left( 0 -\frac{\partial}{\partial x_3} x_2 \right) = 0 \\ r_{33} & = \left( \frac{\partial u_3}{\partial x_3} - \frac{\partial u_3}{\partial x3} \right) = \left( 0 + 0 \right) = 0, \end{align}

and it is calculated as \begin{equation} r_{ij} = \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \end{equation}

Velocity gradient matrix

Finally, we can now calculate the velocity gradient tensor from the calculated strain rate tensor and rotation tensor. This will give as

\begin{align} \frac{\partial u_i}{\partial u_j} = \epsilon_{ij} + \frac{1}{2} r_{ij} = \begin{bmatrix} \frac{1}{1+t} & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{bmatrix} + \frac{1}{2} \begin{bmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix} \end{align}

\begin{equation} \frac{\partial u_i}{\partial u_j} = \begin{bmatrix} \frac{1}{1+t} & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 0 \end{bmatrix}. \end{equation}

Divergence

\begin{align} \nabla \cdot u_{i} & = \partial_i u_i \\ & = \frac{\partial}{x_1} u_1 + \frac{\partial}{x_2} u_2 + \frac{\partial}{x_3} u_3 \\ & = \frac{\partial}{x_1} \frac{x_1}{1+t} + \frac{\partial}{x_2} x_2 + 0 \\ & = \frac{1}{1+t} + 1 \end{align}

\begin{equation} \nabla \cdot u_{i} = \frac{2+t}{1+t} \end{equation}

Vorticity

\begin{align} \nabla \times u_{i} & = \epsilon_{ijk} \partial_j u_k \\ & = \left( \frac{\partial}{\partial x_2} u_3 - \frac{\partial}{\partial x_3} u_2 \right) \hat{x_1} + \left( \frac{\partial}{\partial x_3} u_1 - \frac{\partial}{\partial x_1} u_3 \right) \hat{x_2} + \left( \frac{\partial}{\partial x_1} u_2 - \frac{\partial}{\partial x_2} u_1 \right) \hat{x_3} \\ & = \left( 0 - \frac{\partial}{\partial x_3} x_2 \right) \hat{x_1} + \left( \frac{\partial}{\partial x_3} \frac{x_1}{1+t} - 0 \right) \hat{x_2} + \left( \frac{\partial}{\partial x_1} x_2 - \frac{\partial}{\partial x_2} \frac{x_1}{1+t} \right) \hat{x_3} \\ & = 0 \, \hat{x_1} + 0 \, \hat{x_2} + 0 \, \hat{x_3} \end{align}

\begin{equation} \nabla \times u_{i} = 0 \end{equation}

Simulation of the simple plane flow

def your_u_function(x, y, t):
    u = x / (1+t)
    return u

def your_v_function(x, y, t):
    v = y
    return v

# Define your 2D vector field function (u, v) and time array
def vector_field(x, y, t):
    # Replace with your actual u and v functions
    u = your_u_function(x, y, t)
    v = your_v_function(x, y, t)
    return u, v

# Set up the figure and axes
fig, ax = plt.subplots(figsize=(8, 6))
ax.set_xlim(-2, 2)
ax.set_ylim(-2, 2)

# Create a grid of points
x = np.linspace(-2, 2, 100)
y = np.linspace(-2, 2, 100)
X, Y = np.meshgrid(x, y)

# Function to update the plot for each frame
def update(frame):
    ax.clear()
    ax.set_xlim(-2, 2)
    ax.set_ylim(-2, 2)
    
    # Evaluate the vector field at each point on the grid
    U, V = vector_field(X, Y, frame)
    
    # Plot streamlines
    stream = ax.streamplot(X, Y, U, V, color='black')
    
    # Add a title or other annotations if needed
    ax.set_title(f'Time: {frame}')
    
# Create an animation
num_frames = 15  # Adjust the number of frames as needed
animation = FuncAnimation(fig, update, frames=num_frames, interval=100)

# Save the animation as a GIF
animation.save('streamlines_animation2.gif', writer='pillow')

# Display the animation in the Jupyter Notebook
html_output = HTML(animation.to_jshtml())
display(html_output)

Once Loop Reflect